function out = solver_exploc_matrix(struc_params, solv_params, tau_B, entry_status_exog)
    flag = 0;
    %% get parameters

    % structural parameters
    sig = struc_params.sig; % elasticity of substitution
    alp = struc_params.alp; % input suitability
    m_phi = struc_params.m_phi; % mean of fundamental productivity
    m_del = struc_params.m_del; % mean of fundamental quality
    v_phi = struc_params.v_phi; % std of fundamental productivity
    v_del = struc_params.v_del; % std of fundamental quality
    v_cor = struc_params.v_cor; % correlation between phi and del
    m_psi_II = struc_params.m_psi_II; % scale parameter of psi function
    m_psi_IC = struc_params.m_psi_IC; % scale parameter of psi function
    m_psi_F = struc_params.m_psi_F;
    s_xi = struc_params.s_xi; % shape parameter of xi distribution
    L = struc_params.L; % population size
    tau_D = struc_params.tau_D; % domestic trade cost
    tau_DF = struc_params.tau_DF; % domestic trade cost
    w = struc_params.w; % unit cost of dirty inputs
    t = struc_params.t; % emission tax rate
    eps = struc_params.eps; % scale parameter (carbon content)
    a_LC = struc_params.a_LC; % average relocation costs
    net_open = struc_params.net_open;

    % solver parameters
    qt_min = solv_params.qt_min; % min quantile for phi/del grid
    qt_max = solv_params.qt_max; % max quantile for phi/del grid
    n_grid = solv_params.n_grid; % number of grid points for phi/del grid
    tol = solv_params.tol; % numerical tolerance for iterations
    max_iter = solv_params.max_iter; % max iterations

    mu = sig / (sig - 1); % constant markup
    m_xi = 1; % mean of xi that follows weibull distribution
    lam_xi = m_xi / gamma(1 + 1 / s_xi); % scale parameter of the weibull distribution that xi follows
    G_xi = @(x) (1 - exp(- (x ./ lam_xi) .^ s_xi)) .* (x >= 0); % cdf of xi ~ weibull distribution
    E_xi = @(x) lam_xi * gamma(1 + 1 / s_xi) * ... % partial expectation of xi from 0 to x
        gammainc((x ./ lam_xi) .^ s_xi, 1 + 1 / s_xi);

    n_grid_phi = n_grid;
    n_grid_del = n_grid;

    matching_flag = solv_params.matching_flag;

    if ~strcmp(matching_flag, 'endogenous')

        if isfield(struc_params, 'io_param')

            rng(303, 'twister');
            s = rng;
            rng(s);
            io_param_D = rand(n_grid_phi ^ 2, n_grid_del ^ 2);

            io_param_D = 2 .* io_param_D .* struc_params.io_param_D;

            rng(5280, 'twister');
            sf = rng;
            rng(sf);
            io_param_F = rand(n_grid_phi ^ 2, n_grid_del ^ 2);

            io_param_F = 2 .* io_param_F .* struc_params.io_param_F;

            io_param = struc_params.io_param;
        else
            disp("WARNING NO io_param SUPPLIED USING DEFAULT")
            io_param_D = 0.1232;
            io_param_F = io_param_D * 0.1;
        end

        disp("exogenous solver starting")
    end

    %% Generate random phi and del
    % Fundamental efficiency and quality follow log-normal
    qt_grid = linspace(qt_min, qt_max, n_grid_phi)';
    phi = exp(norminv(qt_grid, m_phi, v_phi));
    out.phi = phi;

    qt_grid = linspace(qt_min, qt_max, n_grid_del)';
    del = exp(norminv(qt_grid, m_del, v_del));
    out.del = del;

    %% Density function g_chi
    M = [m_phi; m_del]; % mean vector of (log(phi),log(del))
    V = [v_phi ^ 2, v_cor * v_phi * v_del; % covariance matrix of (log(phi),log(del))
         v_cor * v_phi * v_del, v_del ^ 2];
    g_chi_fun = @(phi, del) mvnpdf([log(phi); log(del)], M, V) / (phi * del); % density function of (phi,del)

    g_chi = zeros(n_grid_phi, n_grid_del); % density function evaluated at 20x20(phi,del) grid points

    for n_phii = 1:n_grid_phi

        for n_dell = 1:n_grid_del

            if phi(n_phii) == 0 || del(n_dell) == 0
                g_chi(n_phii, n_dell) = 0;
            else
                g_chi(n_phii, n_dell) = g_chi_fun(phi(n_phii), del(n_dell));
            end

        end

    end

    g_chi = g_chi / sum(sum(g_chi)); % normalize the density so that the sum of g_chi is 1

    %% Set the initial values
    % Set initial values of fa and Kappa
    Kappa = ones(n_grid_phi, n_grid_del); % Kappa is the unit cost of utilizing emission-intensive inputs: price of coal + total emission tax / abatement investment
    fa = ones(n_grid_phi, n_grid_del);

    % Set initial values of DelH
    phidelKappa_bar = matmean(((phi) .^ (sig - 1) * ones(1, n_grid_del)) .* (ones(n_grid_phi, 1) * ((del)') .^ (sig - 1)) .* Kappa .^ (-sig), g_chi, 0);
    Del_H_emp = mu ^ sig * L / phidelKappa_bar; % Del_H when empty nework (Del_H is the household demand shifter)
    Del_H_guess = Del_H_emp;

    % Set initial values of Phi and Del
    Phi_guess = zeros(n_grid_phi, n_grid_del);
    Del_guess = zeros(n_grid_phi, n_grid_del);
    phi_matrix = zeros(n_grid_phi, n_grid_del);
    del_matrix = zeros(n_grid_phi, n_grid_del);

    for n = 1:n_grid_del
        Del_guess(:, n) = mu ^ (-sig) * (del(n)) ^ (sig - 1); % Del of inland firms when empty network
        del_matrix(:, n) = del(n);
    end

    out.Del_guess = Del_guess;
    out.del_matrix = del_matrix;

    for n = 1:n_grid_phi
        Phi_guess(n, :) = (phi(n)) ^ (sig - 1);
        phi_matrix(n, :) = phi(n);
    end

    out.Phi_guess = Phi_guess;
    out.phi_matrix = phi_matrix;

    %% Export / Coastal Status Determination
    % Entry status
    entry_status = ones(n_grid_phi, n_grid_del);

    % Coastal / Inland status
    coastal_status = 0.75 .* ones(n_grid_phi, n_grid_del);

    if ~strcmp(matching_flag, 'endogenous')
        coastal_status = 0.75 .* ones(n_grid_phi, n_grid_del);
    end

    inland_status = 1 - coastal_status; % everyone starts inland

    % Total profit of firms after relocation to the coast
    profit_relocation = zeros(n_grid_phi, n_grid_del);

    % Calculate firms' total profits to sort out exiting firms
    prof_shar = 0.0;
    entry_cost = 0.001;
    profit(entry_status, coastal_status, inland_status, 1, 1);
    profit_orig = (1 - prof_shar) .* out.total_profit - entry_cost;
    out.profit_orig = profit_orig;

    % Derive ranking matrix by firms' benchmark (initial) total profits
    A = reshape(profit_orig, 1, []);
    [~, idx] = ismember(A, sort(A, 'descend'));

    FirmRankA = reshape(idx, n_grid_phi, n_grid_del);
    out.FirmRankA = FirmRankA;

    % Derive position matrix by the rank of firms
    idxrow = zeros(n_grid_phi * n_grid_del, 1);
    idxcol = zeros(n_grid_phi * n_grid_del, 1);

    for i = 1:max(idx)
        [idxrow(i, 1), idxcol(i, 1)] = find(FirmRankA == i);
    end

    idxrankA = [idxrow, idxcol];

    out.idxrow = idxrow;
    out.idxcol = idxcol;
    out.idxrank = idxrankA;

    % Calculate total profits of firms after relocating to coast from inland by ranking
    xi_max_profit_relocation = zeros(n_grid_phi, n_grid_del);
    prob_coastal = coastal_status; % zeros(n_grid_phi, n_grid_del);
    prob_coastal(prob_coastal == 1) = 0.99;

    prob_coastal_old = prob_coastal;
    prob_inland = (1 - prob_coastal) .* entry_status;

    % Fix the set of firms to be the same in the exogenous model
    if ~strcmp(matching_flag, 'endogenous')
        entry_status = entry_status_exog;
        prob_coastal(entry_status == 0) = 0;
        prob_inland(entry_status == 0) = 0;
    end

    loc_res = 1000;
    loc_niter = 0;

    while loc_res > tol && loc_niter <= max_iter
        loc_niter = loc_niter + 1;
        disp(loc_niter);

        prob_inland(prob_coastal == 1) = 0.01;
        prob_coastal(prob_coastal == 1) = 0.99;
        prob_coastal(prob_inland == 1) = 0.01;
        prob_inland(prob_inland == 1) = 0.99;

        % RERUN ALL FIRMS LOCATION DECISIONS GIVEN THE PREVIOUS DECISIONS
        % THIS IS THE FINAL RUN
        for rankidx = 1:n_grid_phi * n_grid_del

            if entry_status(idxrow(rankidx), idxcol(rankidx)) == 1

                % firm profit of moving to the coast
                coast_tmp = prob_coastal;
                coast_tmp(idxrow(rankidx), idxcol(rankidx)) = 1;
                inland_tmp = prob_inland;
                inland_tmp(idxrow(rankidx), idxcol(rankidx)) = 0;
                profit_coast = profit(entry_status, coast_tmp, inland_tmp, idxrow(rankidx), idxcol(rankidx));

                % firm profit of remaining inland
                coast_tmp(idxrow(rankidx), idxcol(rankidx)) = 0;
                inland_tmp(idxrow(rankidx), idxcol(rankidx)) = 1;
                profit_inland = profit(entry_status, coast_tmp, inland_tmp, idxrow(rankidx), idxcol(rankidx));

                % Get probablity of locating on the coast for firm-rankidx
                xi_max_profit_relocation(idxrow(rankidx), idxcol(rankidx)) = max((profit_coast - profit_inland) ./ a_LC, 0);

                prob_coastal(idxrow(rankidx), idxcol(rankidx)) = G_xi(xi_max_profit_relocation(idxrow(rankidx), idxcol(rankidx))) .* entry_status(idxrow(rankidx), idxcol(rankidx));
                prob_inland(idxrow(rankidx), idxcol(rankidx)) = (1 - prob_coastal(idxrow(rankidx), idxcol(rankidx))) .* entry_status(idxrow(rankidx), idxcol(rankidx));

                firm_profits = out.total_profit;
                % update the entry status based off firm movement
                profit_tmp = (1 - prof_shar) .* firm_profits - entry_cost;

                if strcmp(matching_flag, 'endogenous')
                    entry_status(profit_tmp < 0) = 0;
                    entry_status(profit_tmp >= 0) = 1;
                end

                prob_coastal(entry_status == 0) = 0;
                prob_inland(entry_status == 0) = 0;
            end

        end

        loc_res = max(max(abs(prob_coastal - prob_coastal_old)));
        % store results for next iteration's residual and continue
        if loc_res > tol
            prob_coastal_old = prob_coastal;
            disp(loc_res);
        end

    end

    if loc_niter == max_iter
        disp("WARNING: RESULTS DID NOT CONVERGE")
    end

    % Calculate the number and share of firms in each region

    N_coastal = matmean(prob_coastal, g_chi, 0);
    N_inland = matmean(prob_inland, g_chi, 0);

    share_coastal = N_coastal / (N_coastal + N_inland);
    share_inland = 1 - share_coastal;

    out.profit_relocation = profit_relocation;
    out.xi_max_profit_relocation = xi_max_profit_relocation;
    out.prob_coastal = prob_coastal;
    out.prob_inland = prob_inland;
    out.N_coastal = N_coastal;
    out.N_inland = N_inland;
    out.share_coastal = share_coastal;
    out.share_inland = share_inland;

    %% Compute and store outcomes
    % [1] Status Calculation
    out.entry_status = entry_status;
    % [2] Network Characteristics by Status
    Phi = out.Phi;
    Del = out.Del;
    Phi_inland = prob_inland .* Phi;
    Phi_coastal = prob_coastal .* Phi;
    Del_inland = prob_inland .* Del;
    Del_coastal = prob_coastal .* Del;

    out.Phi_inland = Phi_inland;
    out.Phi_coastal = Phi_coastal;
    out.Del_inland = Del_inland;
    out.Del_coastal = Del_coastal;

    % [3] Abatement Investment by Status
    fa = entry_status .* fa;
    out.fa = fa;

    total_fa = matmean(fa, g_chi, 0);
    out.total_fa = total_fa;

    fa_inland = entry_status .* prob_inland .* fa;
    fa_coastal = entry_status .* prob_coastal .* fa;

    total_fa_inland = matmean(fa_inland, g_chi, 0);
    total_fa_coastal = matmean(fa_coastal, g_chi, 0);

    out.fa_inland = fa_inland;
    out.fa_coastal = fa_coastal;

    out.total_fa_inland = total_fa_inland;
    out.total_fa_coastal = total_fa_coastal;

    % [3-1] Relative abatement investment compared to revenue
    rev = out.rev;
    rev = entry_status .* rev;

    rel_abate_rev = fa ./ rev;
    rel_abate_rev(isnan(rel_abate_rev)) = 0;
    rel_abate_rev_coastal = rel_abate_rev .* prob_coastal;
    rel_abate_rev_inland = rel_abate_rev .* prob_inland;

    out.rel_abate_rev = rel_abate_rev;
    out.rel_abate_rev_coastal = rel_abate_rev_coastal;
    out.rel_abate_rev_inland = rel_abate_rev_inland;

    avg_rel_abate_rev_coastal = matmean(rel_abate_rev_coastal, g_chi, 0);
    avg_rel_abate_rev_inland = matmean(rel_abate_rev_inland, g_chi, 0);

    out.avg_rel_abate_rev_coastal = avg_rel_abate_rev_coastal;
    out.avg_rel_abate_rev_inland = avg_rel_abate_rev_inland;

    % [4] Unit Cost of Emission-intensive factor by Status
    Kappa = entry_status .* Kappa;
    out.Kappa = Kappa;

    Kappa_inland = entry_status .* prob_inland .* Kappa;
    Kappa_coastal = entry_status .* prob_coastal .* Kappa;
    total_Kappa_inland = matmean(Kappa_inland, g_chi, 0);
    total_Kappa_coastal = matmean(Kappa_coastal, g_chi, 0);

    out.Kappa_inland = Kappa_inland;
    out.Kappa_coastal = Kappa_coastal;
    out.total_Kappa_inland = total_Kappa_inland;
    out.total_Kappa_coastal = total_Kappa_coastal;

    % [5] Key Outputs
    opt_e = out.opt_e;
    opt_int = out.opt_int;
    p_x = out.p_x;
    p_x = entry_status .* p_x;
    X = out.X;

    % Relative Factor Price
    Kappa_px = Kappa ./ p_x;
    px_Kappa = p_x ./ Kappa;

    Kappa_px(isnan(Kappa_px)) = 0;
    px_Kappa(isnan(px_Kappa)) = 0;

    out.Kappa_px = Kappa_px;
    out.px_Kappa = px_Kappa;

    % [6] Total Emissions by Firms
    opt_e_inland = entry_status .* prob_inland .* opt_e;
    opt_e_coastal = entry_status .* prob_coastal .* opt_e;

    out.opt_e_inland = opt_e_inland;
    out.opt_e_coastal = opt_e_coastal;

    % [7] Emission Intensity by Firms
    opt_int_inland = entry_status .* prob_inland .* opt_int;
    opt_int_coastal = entry_status .* prob_coastal .* opt_int;

    out.opt_int_inland = opt_int_inland;
    out.opt_int_coastal = opt_int_coastal;

    % [8] Share of intermediates across inland and coast
    int_x = out.int_x;
    int_x_inland = entry_status .* prob_inland .* int_x;
    int_x_coastal = entry_status .* prob_coastal .* int_x;
    total_int_x_inland = matmean(int_x_inland, g_chi, 0);
    total_int_x_coastal = matmean(int_x_coastal, g_chi, 0);
    share_int_x_inland = total_int_x_inland / (total_int_x_inland + total_int_x_coastal);
    share_int_x_coastal = 1 - share_int_x_inland;

    share_int_x_inland(isnan(share_int_x_inland)) = 0;
    share_int_x_coastal(isnan(share_int_x_coastal)) = 0;

    out.share_int_x_inland = share_int_x_inland;
    out.share_int_x_coastal = share_int_x_coastal;

    % [9] Intermediate to dirty input ratio across inland and coast
    xdratio = out.xdratio;

    xdratio_inland = entry_status .* prob_inland .* xdratio;
    xdratio_coastal = entry_status .* prob_coastal .* xdratio;

    out.xdratio_inland = xdratio_inland;
    out.xdratio_coastal = xdratio_coastal;

    total_xdratio_inland = matmean(xdratio_inland, g_chi, 0);
    total_xdratio_coastal = matmean(xdratio_coastal, g_chi, 0);

    out.total_xdratio_inland = total_xdratio_inland;
    out.total_xdratio_coastal = total_xdratio_coastal;

    dxratio = out.dxratio;
    dxratio_inland = entry_status .* prob_inland .* dxratio;
    dxratio_coastal = entry_status .* prob_coastal .* dxratio;

    out.dxratio_inland = dxratio_inland;
    out.dxratio_coastal = dxratio_coastal;

    % [10] Regional Emission Intensities
    total_X = entry_status .* X;
    inland_X = entry_status .* prob_inland .* X;
    coastal_X = entry_status .* prob_coastal .* X;

    out.total_X = total_X;
    out.inland_X = inland_X;
    out.coastal_X = coastal_X;

    inland_eint = matmean(opt_e_inland, g_chi, 0) / matmean(inland_X, g_chi, 0);
    inland_eint(isnan(inland_eint)) = 0;
    coastal_eint = matmean(opt_e_coastal, g_chi, 0) / matmean(coastal_X, g_chi, 0);
    coastal_eint(isnan(coastal_eint)) = 0;

    out.inland_eint = inland_eint;
    out.coastal_eint = coastal_eint;

    % [11] Total production inland and coast
    total_X_val = matmean(total_X, g_chi, 0);
    total_inland_X = matmean(inland_X, g_chi, 0);
    total_coastal_X = matmean(coastal_X, g_chi, 0);

    out.total_X_val = total_X_val;
    out.total_inland_X = total_inland_X;
    out.total_coastal_X = total_coastal_X;

    % [12] Regional Total Emission
    inland_etotal = matmean(opt_e_inland, g_chi, 0);
    coastal_etotal = matmean(opt_e_coastal, g_chi, 0);

    share_inland_etotal = matmean(opt_e_inland, g_chi, 0) / (matmean(opt_e_inland, g_chi, 0) + matmean(opt_e_coastal, g_chi, 0));
    share_coastal_etotal = 1 - share_inland_etotal;

    out.inland_etotal = inland_etotal;
    out.coastal_etotal = coastal_etotal;

    out.share_inland_etotal = share_inland_etotal;
    out.share_coastal_etotal = share_coastal_etotal;

    % [13] Aggregate emission and emission intensity
    total_emission = matmean(entry_status .* opt_e, g_chi, 0);
    total_eint = matmean(entry_status .* opt_e, g_chi, 0) / matmean(entry_status .* X, g_chi, 0);

    out.total_emission = total_emission;
    out.total_eint = total_eint;

    % [14] Number of International Firms
    prob_intl = (newmatmean(m_C_inland_Finland, g_chi, 0) .* entry_status ...
        + newmatmean(m_C_inland_Fcoastal, g_chi, 0) .* entry_status ...
        + newmatmean(m_C_coastal_Finland, g_chi, 0) .* entry_status ...
        + newmatmean(m_C_coastal_Fcoastal, g_chi, 0) .* entry_status) ...
        + (newmatmean(m_S_inland_Finland, g_chi, 0) .* entry_status ...
        + newmatmean(m_S_inland_Fcoastal, g_chi, 0) .* entry_status ...
        + newmatmean(m_S_coastal_Finland, g_chi, 0) .* entry_status ...
        + newmatmean(m_S_coastal_Fcoastal, g_chi, 0) .* entry_status);

    prob_intl = prob_intl ./ max(max(prob_intl));
    INTL = prob_intl;
    out.INTL = INTL;

    % [15] Number of Domestic Firms
    prob_dom = (1 - prob_intl) .* entry_status;

    DOM = prob_dom;
    out.DOM = DOM;

    % [16] Differentiate firms by their status
    prob_inland_domestic = entry_status .* prob_inland .* DOM;
    prob_inland_exporter = entry_status .* prob_inland .* INTL;
    prob_coastal_domestic = entry_status .* prob_coastal .* DOM;
    prob_coastal_exporter = entry_status .* prob_coastal .* INTL;

    out.prob_inland_domestic = prob_inland_domestic;
    out.prob_inland_exporter = prob_inland_exporter;
    out.prob_coastal_domestic = prob_coastal_domestic;
    out.prob_coastal_exporter = prob_coastal_exporter;

    N_inland_domestic = matmean(prob_inland_domestic, g_chi, 0);
    N_inland_exporter = matmean(prob_inland_exporter, g_chi, 0);
    N_coastal_domestic = matmean(prob_coastal_domestic, g_chi, 0);
    N_coastal_exporter = matmean(prob_coastal_exporter, g_chi, 0);

    out.N_inland_domestic = N_inland_domestic;
    out.N_inland_exporter = N_inland_exporter;
    out.N_coastal_domestic = N_coastal_domestic;
    out.N_coastal_exporter = N_coastal_exporter;

    share_inland_domestic = N_inland_domestic / (N_inland_domestic + N_inland_exporter + N_coastal_domestic + N_coastal_exporter);
    share_inland_exporter = N_inland_exporter / (N_inland_domestic + N_inland_exporter + N_coastal_domestic + N_coastal_exporter);
    share_coastal_domestic = N_coastal_domestic / (N_inland_domestic + N_inland_exporter + N_coastal_domestic + N_coastal_exporter);
    share_coastal_exporter = N_coastal_exporter / (N_inland_domestic + N_inland_exporter + N_coastal_domestic + N_coastal_exporter);

    out.share_inland_domestic = share_inland_domestic;
    out.share_inland_exporter = share_inland_exporter;
    out.share_coastal_domestic = share_coastal_domestic;
    out.share_coastal_exporter = share_coastal_exporter;

    share_dom = (N_inland_domestic + N_coastal_domestic) / (N_inland_domestic + N_inland_exporter + N_coastal_domestic + N_coastal_exporter);
    share_intl = (N_inland_exporter + N_coastal_exporter) / (N_inland_domestic + N_inland_exporter + N_coastal_domestic + N_coastal_exporter);

    out.share_dom = share_dom;
    out.share_intl = share_intl;

    % [17] Relative price of inputs
    relprice_dx = Kappa ./ p_x;
    relprice_xd = p_x ./ Kappa;
    relprice_dx(isnan(relprice_dx)) = 0;
    relprice_xd(isnan(relprice_xd)) = 0;

    relprice_dx_inland = entry_status .* prob_inland .* relprice_dx;
    relprice_dx_coastal = entry_status .* prob_coastal .* relprice_dx;
    relprice_xd_inland = entry_status .* prob_inland .* relprice_xd;
    relprice_xd_coastal = entry_status .* prob_coastal .* relprice_xd;

    out.relprice_dx = relprice_dx;
    out.relprice_xd = relprice_xd;

    out.relprice_dx_inland = relprice_dx_inland;
    out.relprice_dx_coastal = relprice_dx_coastal;
    out.relprice_xd_inland = relprice_xd_inland;
    out.relprice_xd_coastal = relprice_xd_coastal;

    % [18] Abatment investment vs. Outsourcing
    abate_outsource_ratio = fa ./ (int_x);
    outsource_abate_ratio = (int_x) ./ fa;

    abate_outsource_ratio(isnan(abate_outsource_ratio)) = 0;
    outsource_abate_ratio(isnan(outsource_abate_ratio)) = 0;

    out.abate_outsource_ratio = abate_outsource_ratio;
    out.outsource_abate_ratio = outsource_abate_ratio;

    abate_outsource_ratio_inland = entry_status .* prob_inland .* abate_outsource_ratio;
    abate_outsource_ratio_coastal = entry_status .* prob_coastal .* abate_outsource_ratio;

    outsource_abate_ratio_inland = entry_status .* prob_inland .* outsource_abate_ratio;
    outsource_abate_ratio_coastal = entry_status .* prob_coastal .* outsource_abate_ratio;

    out.abate_outsource_ratio_inland = abate_outsource_ratio_inland;
    out.abate_outsource_ratio_coastal = abate_outsource_ratio_coastal;
    out.outsource_abate_ratio_inland = outsource_abate_ratio_inland;
    out.outsource_abate_ratio_coastal = outsource_abate_ratio_coastal;

    % [19] Abatment investment vs. Outsourcing
    abate_outsourceexp_ratio = fa ./ (int_x .* p_x);
    outsourceexp_abate_ratio = (int_x .* p_x) ./ fa;

    abate_outsourceexp_ratio(isnan(abate_outsourceexp_ratio)) = 0;
    outsourceexp_abate_ratio(isnan(outsourceexp_abate_ratio)) = 0;

    out.abate_outsourceexp_ratio = abate_outsourceexp_ratio;
    out.outsourceexp_abate_ratio = outsourceexp_abate_ratio;

    abate_outsourceexp_ratio_inland = entry_status .* prob_inland .* abate_outsourceexp_ratio;
    abate_outsourceexp_ratio_coastal = entry_status .* prob_coastal .* abate_outsourceexp_ratio;

    outsourceexp_abate_ratio_inland = entry_status .* prob_inland .* outsourceexp_abate_ratio;
    outsourceexp_abate_ratio_coastal = entry_status .* prob_coastal .* outsourceexp_abate_ratio;

    out.abate_outsourceexp_ratio_inland = abate_outsourceexp_ratio_inland;
    out.abate_outsourceexp_ratio_coastal = abate_outsourceexp_ratio_coastal;
    out.outsourceexp_abate_ratio_inland = outsourceexp_abate_ratio_inland;
    out.outsourceexp_abate_ratio_coastal = outsourceexp_abate_ratio_coastal;

    % [20] Share of abatement investment from total exp to reduce emissions
    abate_totalexp_ratio = fa ./ (fa + int_x .* p_x);
    outsourceexp_totalexp_ratio = (int_x .* p_x) ./ (fa + int_x .* p_x);
    abate_totalexp_ratio(isnan(abate_totalexp_ratio)) = 0;
    outsourceexp_totalexp_ratio(isnan(outsourceexp_totalexp_ratio)) = 0;

    out.abate_totalexp_ratio = abate_totalexp_ratio;
    out.outsourceexp_totalexp_ratio = outsourceexp_totalexp_ratio;

    abate_totalexp_ratio_inland = entry_status .* prob_inland .* abate_totalexp_ratio;
    abate_totalexp_ratio_coastal = entry_status .* prob_coastal .* abate_totalexp_ratio;

    outsourceexp_totalexp_ratio_inland = entry_status .* prob_inland .* outsourceexp_totalexp_ratio;
    outsourceexp_totalexp_ratio_coastal = entry_status .* prob_coastal .* outsourceexp_totalexp_ratio;

    out.abate_totalexp_ratio_inland = abate_totalexp_ratio_inland;
    out.abate_totalexp_ratio_coastal = abate_totalexp_ratio_coastal;
    out.outsourceexp_totalexp_ratio_inland = outsourceexp_totalexp_ratio_inland;
    out.outsourceexp_totalexp_ratio_coastal = outsourceexp_totalexp_ratio_coastal;

    % [21] Input use by firm types
    d = out.d;
    d_coastal = entry_status .* prob_coastal .* d;
    d_inland = entry_status .* prob_inland .* d;
    d_coastal_domestic = entry_status .* prob_coastal_domestic .* d;
    d_coastal_exporter = entry_status .* prob_coastal_exporter .* d;
    d_inland_domestic = entry_status .* prob_inland_domestic .* d;
    d_inland_exporter = entry_status .* prob_inland_exporter .* d;

    out.d_coastal = d_coastal;
    out.d_inland = d_inland;
    out.d_coastal_domestic = d_coastal_domestic;
    out.d_coastal_exporter = d_coastal_exporter;
    out.d_inland_domestic = d_inland_domestic;
    out.d_inland_exporter = d_inland_exporter;

    int_x = out.int_x;
    int_x_coastal = entry_status .* prob_coastal .* int_x;
    int_x_inland = entry_status .* prob_inland .* int_x;
    int_x_coastal_domestic = entry_status .* prob_coastal_domestic .* int_x;
    int_x_coastal_exporter = entry_status .* prob_coastal_exporter .* int_x;
    int_x_inland_domestic = entry_status .* prob_inland_domestic .* int_x;
    int_x_inland_exporter = entry_status .* prob_inland_exporter .* int_x;

    out.int_x_coastal = int_x_coastal;
    out.int_x_inland = int_x_inland;
    out.int_x_coastal_domestic = int_x_coastal_domestic;
    out.int_x_coastal_exporter = int_x_coastal_exporter;
    out.int_x_inland_domestic = int_x_inland_domestic;
    out.int_x_inland_exporter = int_x_inland_exporter;

    % [22] d/fa ratio
    d_fa = entry_status .* (d ./ fa);
    d_fa(isnan(d_fa)) = 0;

    out.d_fa = d_fa;
    d_fa_coastal = entry_status .* prob_coastal .* d_fa;
    d_fa_inland = entry_status .* prob_inland .* d_fa;
    d_fa_coastal_domestic = entry_status .* prob_coastal_domestic .* d_fa;
    d_fa_caostal_exporter = entry_status .* prob_coastal_exporter .* d_fa;
    d_fa_inland_domestic = entry_status .* prob_inland_domestic .* d_fa;
    d_fa_inland_exporter = entry_status .* prob_inland_exporter .* d_fa;

    out.d_fa_coastal = d_fa_coastal;
    out.d_fa_inland = d_fa_inland;
    out.d_fa_coastal_domestic = d_fa_coastal_domestic;
    out.d_fa_caostal_exporter = d_fa_caostal_exporter;
    out.d_fa_inland_domestic = d_fa_inland_domestic;
    out.d_fa_inland_exporter = d_fa_inland_exporter;

    % [23] fa/d ratio
    fa_d = entry_status .* (fa ./ d);
    fa_d(isnan(fa_d)) = 0;

    out.fa_d = fa_d;

    fa_d_coastal = entry_status .* prob_coastal .* fa_d;
    fa_d_inland = entry_status .* prob_inland .* fa_d;
    fa_d_coastal_domestic = entry_status .* prob_coastal_domestic .* fa_d;
    fa_d_caostal_exporter = entry_status .* prob_coastal_exporter .* fa_d;
    fa_d_inland_domestic = entry_status .* prob_inland_domestic .* fa_d;
    fa_d_inland_exporter = entry_status .* prob_inland_exporter .* fa_d;

    out.fa_d_coastal = fa_d_coastal;
    out.fa_d_inland = fa_d_inland;
    out.fa_d_coastal_domestic = fa_d_coastal_domestic;
    out.fa_d_caostal_exporter = fa_d_caostal_exporter;
    out.fa_d_inland_domestic = fa_d_inland_domestic;
    out.fa_d_inland_exporter = fa_d_inland_exporter;

    % [24] int_x share out of total input use
    share_intx = (int_x ./ (d + int_x));
    share_intx(isnan(share_intx)) = 0;
    share_intx = entry_status .* share_intx;

    share_intx_inland = prob_inland .* share_intx;
    share_intx_coastal = prob_coastal .* share_intx;

    share_intx_inland_domestic = prob_inland_domestic .* share_intx;
    share_intx_coastal_domestic = prob_coastal_domestic .* share_intx;
    share_intx_inland_exporter = prob_inland_exporter .* share_intx;
    share_intx_coastal_exporter = prob_coastal_exporter .* share_intx;

    out.share_intx = share_intx;
    out.share_intx_inland = share_intx_inland;
    out.share_intx_coastal = share_intx_coastal;
    out.share_intx_inland_domestic = share_intx_inland_domestic;
    out.share_intx_coastal_domestic = share_intx_coastal_domestic;
    out.share_intx_inland_exporter = share_intx_inland_exporter;
    out.share_intx_coastal_exporter = share_intx_coastal_exporter;

    % [25]  d share out of total input use
    share_d = (d ./ (d + int_x));
    share_d(isnan(share_d)) = 0;
    share_d = entry_status .* share_d;

    share_d_inland = prob_inland .* share_d;
    share_d_coastal = prob_coastal .* share_d;

    share_d_inland_domestic = prob_inland_domestic .* share_d;
    share_d_coastal_domestic = prob_coastal_domestic .* share_d;
    share_d_inland_exporter = prob_inland_exporter .* share_d;
    share_d_coastal_exporter = prob_coastal_exporter .* share_d;

    out.share_d = share_d;
    out.share_d_inland = share_d_inland;
    out.share_d_coastal = share_d_coastal;
    out.share_d_inland_domestic = share_d_inland_domestic;
    out.share_d_coastal_domestic = share_d_coastal_domestic;
    out.share_d_inland_exporter = share_d_inland_exporter;
    out.share_d_coastal_exporter = share_d_coastal_exporter;

    % [26] int_x expenditure share out of total expenditure on reducing emissions
    share_intxexp = ((int_x .* p_x) ./ (fa + int_x .* p_x));
    share_intxexp(isnan(share_intxexp)) = 0;
    share_intxexp = entry_status .* share_intxexp;

    share_intxexp_inland = prob_inland .* share_intxexp;
    share_intxexp_coastal = prob_coastal .* share_intxexp;

    share_intxexp_inland_domestic = prob_inland_domestic .* share_intxexp;
    share_intxexp_coastal_domestic = prob_coastal_domestic .* share_intxexp;
    share_intxexp_inland_exporter = prob_inland_exporter .* share_intxexp;
    share_intxexp_coastal_exporter = prob_coastal_exporter .* share_intxexp;

    out.share_intxexp = share_intxexp;
    out.share_intxexp_inland = share_intxexp_inland;
    out.share_intxexp_coastal = share_intxexp_coastal;
    out.share_intxexp_inland_domestic = share_intxexp_inland_domestic;
    out.share_intxexp_coastal_domestic = share_intxexp_coastal_domestic;
    out.share_intxexp_inland_exporter = share_intxexp_inland_exporter;
    out.share_intxexp_coastal_exporter = share_intxexp_coastal_exporter;

    % [27] fa expenditure share out of total expenditure on reducing emissions
    share_fa = ((fa) ./ (fa + int_x .* p_x));
    share_fa(isnan(share_fa)) = 0;
    share_fa = entry_status .* share_fa;

    share_fa_inland = prob_inland .* share_fa;
    share_fa_coastal = prob_coastal .* share_fa;

    share_fa_inland_domestic = prob_inland_domestic .* share_fa;
    share_fa_coastal_domestic = prob_coastal_domestic .* share_fa;
    share_fa_inland_exporter = prob_inland_exporter .* share_fa;
    share_fa_coastal_exporter = prob_coastal_exporter .* share_fa;

    out.share_fa = share_fa;
    out.share_fa_inland = share_fa_inland;
    out.share_fa_coastal = share_fa_coastal;
    out.share_fa_inland_domestic = share_fa_inland_domestic;
    out.share_fa_coastal_domestic = share_fa_coastal_domestic;
    out.share_fa_inland_exporter = share_fa_inland_exporter;
    out.share_fa_coastal_exporter = share_fa_coastal_exporter;

    % [28]  MATCHING PROBS
    if matching_flag == "endogenous"
        m_S_inland_inland = matmean(newmatmean(m_S_inland_inland, g_chi, 0), g_chi, 0);
        m_S_inland_coastal = matmean(newmatmean(m_S_inland_coastal, g_chi, 0), g_chi, 0);
        m_S_inland_Finland = matmean(newmatmean(m_S_inland_Finland, g_chi, 0), g_chi, 0);
        m_S_inland_Fcoastal = matmean(newmatmean(m_S_inland_Fcoastal, g_chi, 0), g_chi, 0);
        m_S_coastal_inland = matmean(newmatmean(m_S_coastal_inland, g_chi, 0), g_chi, 0);
        m_S_coastal_coastal = matmean(newmatmean(m_S_coastal_coastal, g_chi, 0), g_chi, 0);
        m_S_coastal_Finland = matmean(newmatmean(m_S_coastal_Finland, g_chi, 0), g_chi, 0);
        m_S_coastal_Fcoastal = matmean(newmatmean(m_S_coastal_Fcoastal, g_chi, 0), g_chi, 0);

        m_C_inland_inland = matmean(newmatmean(m_C_inland_inland, g_chi, 0), g_chi, 0);
        m_C_inland_coastal = matmean(newmatmean(m_C_inland_coastal, g_chi, 0), g_chi, 0);
        m_C_inland_Finland = matmean(newmatmean(m_C_inland_Finland, g_chi, 0), g_chi, 0);
        m_C_inland_Fcoastal = matmean(newmatmean(m_C_inland_Fcoastal, g_chi, 0), g_chi, 0);
        m_C_coastal_inland = matmean(newmatmean(m_C_coastal_inland, g_chi, 0), g_chi, 0);
        m_C_coastal_coastal = matmean(newmatmean(m_C_coastal_coastal, g_chi, 0), g_chi, 0);
        m_C_coastal_Finland = matmean(newmatmean(m_C_coastal_Finland, g_chi, 0), g_chi, 0);
        m_C_coastal_Fcoastal = matmean(newmatmean(m_C_coastal_Fcoastal, g_chi, 0), g_chi, 0);
    else
        m_S_inland_inland = io_param_D;
        m_S_inland_inland(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
        m_S_inland_coastal = io_param_D;
        m_S_inland_coastal(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
        m_S_inland_Finland = io_param_F;
        m_S_inland_Finland(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
        m_S_inland_Fcoastal = io_param_F;
        m_S_inland_Fcoastal(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
        m_S_coastal_inland = io_param_D;
        m_S_coastal_inland(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
        m_S_coastal_coastal = io_param_D;
        m_S_coastal_coastal(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
        m_S_coastal_Finland = io_param_F;
        m_S_coastal_Finland(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
        m_S_coastal_Fcoastal = io_param_F;
        m_S_coastal_Fcoastal(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;

        m_C_inland_inland = io_param_D;
        m_C_inland_inland(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
        m_C_inland_coastal = io_param_D;
        m_C_inland_coastal(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
        m_C_inland_Finland = io_param_F;
        m_C_inland_Finland(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
        m_C_inland_Fcoastal = io_param_F;
        m_C_inland_Fcoastal(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
        m_C_coastal_inland = io_param_D;
        m_C_coastal_inland(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
        m_C_coastal_coastal = io_param_D;
        m_C_coastal_coastal(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
        m_C_coastal_Finland = io_param_F;
        m_C_coastal_Finland(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
        m_C_coastal_Fcoastal = io_param_F;
        m_C_coastal_Fcoastal(kron(entry_status, ones(n_grid_phi, n_grid_del)) == 0) = 0;
    end

    avg_m_S_d = matmean((newmatmean(m_S_inland_inland, g_chi, 0) ...
        + newmatmean(m_S_inland_coastal, g_chi, 0)) ...
        + (newmatmean(m_S_coastal_inland, g_chi, 0) ...
        + newmatmean(m_S_coastal_coastal, g_chi, 0)), g_chi, 0);

    avg_m_C_d = matmean((newmatmean(m_C_inland_inland, g_chi, 0) ...
        + newmatmean(m_C_inland_coastal, g_chi, 0)) ...
        + (newmatmean(m_C_coastal_inland, g_chi, 0) ...
        + newmatmean(m_C_coastal_coastal, g_chi, 0)), g_chi, 0);

    avg_m_C_F = matmean((newmatmean(m_C_inland_Finland, g_chi, 0) ...
        + newmatmean(m_C_inland_Fcoastal, g_chi, 0)) ...
        + (newmatmean(m_C_coastal_Finland, g_chi, 0) ...
        + newmatmean(m_C_coastal_Fcoastal, g_chi, 0)), g_chi, 0);

    avg_m_S_F = matmean((newmatmean(m_S_inland_Finland, g_chi, 0) ...
        + newmatmean(m_S_inland_Fcoastal, g_chi, 0)) ...
        + (newmatmean(m_S_coastal_Finland, g_chi, 0) ...
        + newmatmean(m_S_coastal_Fcoastal, g_chi, 0)), g_chi, 0);

    out.m_S_inland_inland = m_S_inland_inland;
    out.m_S_inland_coastal = m_S_inland_coastal;
    out.m_S_coastal_inland = m_S_coastal_inland;
    out.m_S_coastal_coastal = m_S_coastal_coastal;

    out.m_S_inland_Finland = m_S_inland_Finland;
    out.m_S_inland_Fcoastal = m_S_inland_Fcoastal;
    out.m_S_coastal_Finland = m_S_coastal_Finland;
    out.m_S_coastal_Fcoastal = m_S_coastal_Fcoastal;

    out.m_C_inland_inland = m_C_inland_inland;
    out.m_C_inland_coastal = m_C_inland_coastal;
    out.m_C_coastal_inland = m_C_coastal_inland;
    out.m_C_coastal_coastal = m_C_coastal_coastal;

    out.m_C_inland_Finland = m_C_inland_Finland;
    out.m_C_inland_Fcoastal = m_C_inland_Fcoastal;
    out.m_C_coastal_Finland = m_C_coastal_Finland;
    out.m_C_coastal_Fcoastal = m_C_coastal_Fcoastal;

    mean_matchingprob_inland = (matmean(newmatmean(m_S_inland_inland, g_chi, 0), g_chi, 0) + ...
        matmean(newmatmean(m_S_inland_Finland, g_chi, 0), g_chi, 0) + ...
        matmean(newmatmean(m_S_inland_coastal, g_chi, 0), g_chi, 0) + ...
        matmean(newmatmean(m_S_inland_Fcoastal, g_chi, 0), g_chi, 0)); % .* ones(n_grid_phi, n_grid_del);
    mean_matchingprob_coastal = (matmean(newmatmean(m_S_coastal_inland, g_chi, 0), g_chi, 0) + ...
        matmean(newmatmean(m_S_coastal_Finland, g_chi, 0), g_chi, 0) + ...
        matmean(newmatmean(m_S_coastal_coastal, g_chi, 0), g_chi, 0) + ...
        matmean(newmatmean(m_S_coastal_Fcoastal, g_chi, 0), g_chi, 0));

    mean_matchingprob = mean([mean_matchingprob_coastal, mean_matchingprob_inland]);

    if matching_flag == "exogenous"
        mean_matchingprob = io_param;
    end

    out.mean_matchingprob = mean_matchingprob; % all elements are 0.1265
    out.avg_m_S_d = avg_m_S_d;
    out.avg_m_C_d = avg_m_C_d;
    out.avg_m_S_F = avg_m_S_F;
    out.avg_m_C_F = avg_m_C_F;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Subfunctions %%

    %[1] Definition of profit function
    function total_profit = profit(entry_status, coastal_status, inland_status, i, j)
        % Del_H, PhiDel, fa loop Start Here %%%%%%%%%%%%%%%%
        Del_H = Del_H_guess;
        Del_H_niter = 0;
        Del_H_res = Inf;
        % Del_H loop locates outer than PhiDel loop b/c it is determined after Phi (firms' productivity from network) and Del (demand of firms' product: network demand) is chosen
        while Del_H_res > tol && Del_H_niter <= max_iter
            % increment and display iteration number
            Del_H_niter = Del_H_niter + 1;
            % initialize guess of Phi and Del on first iteration
            if Del_H_niter == 1
                Phi = Phi_guess;
                Del = Del_guess;
            end

            out.Phi_guess = Phi_guess;
            out.Del_guess = Del_guess;

            %% PhiDel Loop Begins
            PhiDel_res = Inf;
            PhiDel_iter = 0;

            while PhiDel_res > tol && PhiDel_iter <= max_iter
                % Increment iteration counter

                PhiDel_iter = PhiDel_iter + 1;
                % iterate on (Phi,Del) functional equations to compute new guess of
                % Phi and Del

                out_PhiDel = PhiDelIter(entry_status, coastal_status, inland_status, Phi, Del, Del_H);

                Phi_inland = out_PhiDel.Phi_inland;
                Del_inland = out_PhiDel.Del_inland;
                Phi_coastal = out_PhiDel.Phi_coastal;
                Del_coastal = out_PhiDel.Del_coastal;

                Phi_new_inland = out_PhiDel.Phi_new_inland;
                Del_new_inland = out_PhiDel.Del_new_inland;
                Phi_new_coastal = out_PhiDel.Phi_new_coastal;
                Del_new_coastal = out_PhiDel.Del_new_coastal;

                RC_inland = out_PhiDel.RC_inland;
                RC_coastal = out_PhiDel.RC_coastal;
                RC_total = out_PhiDel.RC_total;

                % compute PhiDel residual

                Phi_res_inland = max(max(abs(Phi_inland - Phi_new_inland)));
                Del_res_inland = max(max(abs(Del_inland - Del_new_inland)));
                Phi_res_coastal = max(max(abs(Phi_coastal - Phi_new_coastal)));
                Del_res_coastal = max(max(abs(Del_coastal - Del_new_coastal)));

                PhiDel_res_inland = max(Phi_res_inland, Del_res_inland);
                PhiDel_res_coastal = max(Phi_res_coastal, Del_res_coastal);

                PhiDel_res = max(PhiDel_res_inland, PhiDel_res_coastal);

                % update guess of Phi and Del if needed
                if PhiDel_res > tol
                    Phi = Phi_new_inland + Phi_new_coastal;
                    Del = Del_new_inland + Del_new_coastal;

                end

                % Start fa loop here
                fa_res = Inf;
                fa_iter = 0;

                while fa_res > tol && fa_iter <= max_iter \
                    fa_iter = fa_iter + 1;
                    % iterate on fa functional equations to compute new guess of
                    % fa
                    out_fa = faIter(fa, Del, Del_H);
                    fa_new = out_fa.fa_new;
                    % compute fa residual
                    fa_res = max(max(abs(fa - fa_new)));
                    % update guess of fa if needed
                    if fa_res > tol
                        fa = fa_new;
                    end
                end % End of fa loop
                % Calculate the marginal cost of emission-intensive input
                Kappa = w + ((t * eps) ./ (fa));
                Kappa(Kappa > 10) = 4.5; % Cap the max in order for plotting graphs
                Kappa(Kappa == 0) = NaN;
            end % End of PhiDel loop

            Del_H_new = max(L / (matmean(Del .* Kappa .^ (-sig) .* (phi .^ (sig - 1) * ones(1, n_grid_del)), g_chi, 0)), 0);

            % compute and display Del_H_residual
            Del_H_res = abs(Del_H - Del_H_new);

            % update guess of Del_H
            if Del_H_res > tol
                Del_H = Del_H_new;
            end
        end % End of DelH loop

        out.Del_H = Del_H;
        out.Phi = Phi;
        out.Del = Del;
        out.RC_total = RC_total;
        out.fa = fa;
        out.Kappa = Kappa;
        out.phi = phi;
        out.del = del;

        Del_inland = entry_status .* inland_status .* Del;
        Phi_inland = entry_status .* inland_status .* Phi;
        Del_coastal = entry_status .* coastal_status .* Del;
        Phi_coastal = entry_status .* coastal_status .* Phi;

        Pi_inland = (mu - 1) * Del_H .* Del_inland .* Phi_inland - RC_inland;
        Pi_coastal = (mu - 1) * Del_H .* Del_coastal .* Phi_coastal - RC_coastal;

        % Calculate firms' total revenue
        rev = mu * Del_H * Phi .* Del;

        % Calculate firms' relocation costs
        reloc_cost = a_LC * rev;

        % Derive firms' total profits
        out.total_profit = entry_status .* inland_status .* Pi_inland ...
            + entry_status .* coastal_status .* Pi_coastal;

        total_profit = entry_status(i, j) * inland_status(i, j) * Pi_inland(i, j) ...
            + entry_status(i, j) * coastal_status(i, j) * Pi_coastal(i, j);

        % store outputs
        out.reloc_cost = reloc_cost;
        out.rev = rev;

        % matrix of variable profits by firms
        var_prof = (mu - 1) * Del_H * Del .* Phi;
        out.var_prof = var_prof;

        % calculate the optimal emissions
        opt_e = eps * (phi .^ (sig - 1) * ones(1, n_grid_del)) .* Del_H .* Del .* ((fa * w + t * eps) .* fa .^ ((1 - sig) / sig)) .^ (-sig);

        % calculate the optimal firm emission intensity
        PHI = Phi .^ (-sig / (sig - 1));
        PHI(~isfinite(PHI)) = 0; % replace Inf to 0; Inf is from 0 (exit) in Phi
        opt_int = eps * (phi .^ (sig - 1) * ones(1, n_grid_del)) .* PHI .* ((fa * w + t * eps) .* fa .^ ((1 - sig) / sig)) .^ (-sig);
        out.opt_e = opt_e;
        out.opt_int = opt_int;

        % Price index and Utility
        % CPI
        P_H = mu * matmean(Phi .* (ones(n_grid_phi, 1) * ((del) .^ (sig - 1))'), g_chi, 0) ^ (1 / (1 - sig));

        % welfare
        U = Del_H * P_H ^ (-sig);

        % firm's total output production
        X = Del_H * Del .* Phi ^ (sig / (sig - 1));
        X = real(X); % drop imaginary numbers which are very close to zero

        % emission-intensive factor usage
        d = Del_H * Del .* Kappa .^ (-sig) .* kron((phi) .^ (sig - 1), ones(1, n_grid_del));
        total_d = matmean(d, g_chi, 0);
        out.total_d = total_d;

        %%% load matching probabilities %%%
        m_S_inland_inland = out_PhiDel.m_S_inland_inland;
        m_S_inland_coastal = out_PhiDel.m_S_inland_coastal;
        m_S_coastal_inland = out_PhiDel.m_S_coastal_inland;
        m_S_coastal_coastal = out_PhiDel.m_S_coastal_coastal;
        m_S_inland_Finland = out_PhiDel.m_S_inland_Finland;
        m_S_inland_Fcoastal = out_PhiDel.m_S_inland_Fcoastal;
        m_S_coastal_Finland = out_PhiDel.m_S_coastal_Finland;
        m_S_coastal_Fcoastal = out_PhiDel.m_S_coastal_Fcoastal;

        m_S_d1 = m_S_inland_inland;
        m_S_d2 = m_S_inland_coastal;
        m_S_d3 = m_S_coastal_inland;
        m_S_d4 = m_S_coastal_coastal;
        m_S_d = m_S_d1 + m_S_d2 +m_S_d3 + m_S_d4;

        m_S_f1 = m_S_inland_Finland;
        m_S_f2 = m_S_inland_Fcoastal;
        m_S_f3 = m_S_coastal_Finland;
        m_S_f4 = m_S_coastal_Fcoastal;
        m_S_f = m_S_f1 + m_S_f2 +m_S_f3 + m_S_f4;

        out.m_S_d1 = m_S_d1;
        out.m_S_d2 = m_S_d2;
        out.m_S_d3 = m_S_d3;
        out.m_S_d4 = m_S_d4;
        out.m_S_d = m_S_d;
        out.m_S_f1 = m_S_f1;
        out.m_S_f2 = m_S_f2;
        out.m_S_f3 = m_S_f3;
        out.m_S_f4 = m_S_f4;
        out.m_S_f = m_S_f;

        % Find each firms' avg. matching
        m = newmatmean(m_S_d, g_chi, 0) + newmatmean(m_S_f, g_chi, 0);
        p_x = mu .* (Phi) .^ (1 / (1 - sig));
        p_x(~isfinite(p_x)) = 0;
        int_x = alp ^ (sig - 1) * mu ^ (-sig) * Del_H * Del .* (Phi) .^ (sig / (sig - 1));
        xratio = (int_x .* p_x) ./ ((int_x .* p_x) + (d .* Kappa));
        xratio(~isfinite(xratio)) = 0;

        % relative cost share of dirty inputs
        dratio = (d .* Kappa) ./ ((int_x .* p_x) + (d .* Kappa));
        dratio(isnan(dratio)) = 0;

        % factor intensity of firms
        xdratio = int_x ./ d;
        xdratio(isnan(xdratio)) = 0;
        out.xdratio = xdratio;
        dxratio = d ./ int_x;
        dxratio(isnan(dxratio)) = 0;
        out.dxratio = dxratio;

        % store outputs
        out.P_H = P_H;
        out.U = U;
        out.X = X;
        out.d = d;
        out.int_x = int_x;
        out.p_x = p_x;
        out.xratio = xratio;
        out.dratio = dratio;
        out.m = m;

    end % profit function ends

    % [2] Definition of PhiDelIter function
    function out = PhiDelIter(entry_status_iter, coastal_status_iter, inland_status_iter, Phi_iter, Del_iter, Del_H_iter)
        if flag == 1
            disp(max(max(entry_status_iter)));
        end

        Del_inland = entry_status_iter .* inland_status_iter .* Del_iter;
        Phi_inland = entry_status_iter .* inland_status_iter .* Phi_iter;

        Del_coastal = entry_status_iter .* coastal_status_iter .* Del_iter;
        Phi_coastal = entry_status_iter .* coastal_status_iter .* Phi_iter;

        % a matrix of bilateral profits of inland firms from matching with all supply firms in inland/coastal region of home country (n_grid_phi^2 x n_grid_del^2)
        pi_S_inland_inland = (net_open) .* mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Del_inland, Phi_inland);
        pi_S_inland_coastal = (net_open) .* (tau_D) ^ (1 - sig) * mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Del_inland, Phi_coastal);

        % a matrix of bilateral profits of inland firms from matching with all supply firms in inland/coastal region of foreign country
        pi_S_inland_Finland = (net_open) .* (tau_B * tau_D * tau_DF) ^ (1 - sig) * mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Del_inland, Phi_inland);
        pi_S_inland_Fcoastal = (net_open) .* (tau_B * tau_D) ^ (1 - sig) * mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Del_inland, Phi_coastal);

        % a matrix of bilateral profits of coastal firms from matching with all supply firms in inland/coastal region of home country
        pi_S_coastal_inland = (net_open) .* (tau_D) ^ (1 - sig) * mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Del_coastal, Phi_inland);
        pi_S_coastal_coastal = (net_open) .* mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Del_coastal, Phi_coastal);

        % a matrix of bilateral profits of coastal firms from matching with all supply firms in inland/coastal region of foreign country
        pi_S_coastal_Finland = (net_open) .* (tau_B * tau_DF) ^ (1 - sig) * mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Del_coastal, Phi_inland);
        pi_S_coastal_Fcoastal = (net_open) .* (tau_B) ^ (1 - sig) * mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Del_coastal, Phi_coastal);

        % a matrix of bilateral profits of inland firms from matching with all customer firms in inland/coastal region of home country
        pi_C_inland_inland = (net_open) .* mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Phi_inland, Del_inland);
        pi_C_inland_coastal = (net_open) .* (tau_D) ^ (1 - sig) * mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Phi_inland, Del_coastal);

        % a matrix of bilateral profits of inland firms from matching with all customer firms in inland/coastal region of foreign country
        pi_C_inland_Finland = (net_open) .* (tau_B * tau_D * tau_DF) ^ (1 - sig) * mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Phi_inland, Del_inland);
        pi_C_inland_Fcoastal = (net_open) .* (tau_B * tau_D) ^ (1 - sig) * mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Phi_inland, Del_coastal);

        % a matrix of bilateral profits of coastal firms from matching with all customer firms in inland/coastal region of home country
        pi_C_coastal_inland = (net_open) .* (tau_D) ^ (1 - sig) * mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Phi_coastal, Del_inland);
        pi_C_coastal_coastal = (net_open) .* mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Phi_coastal, Del_coastal);

        % a matrix of bilateral profits of coastal firms from matching with all customer firms in inland/coastal region of foreign country
        pi_C_coastal_Finland = (net_open) .* (tau_B * tau_DF) ^ (1 - sig) * mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Phi_coastal, Del_inland);
        pi_C_coastal_Fcoastal = (net_open) .* (tau_B) ^ (1 - sig) * mu ^ (-sig) * (mu - 1) * alp ^ (sig - 1) * Del_H_iter * kron(Phi_coastal, Del_coastal);

        % Average relationship costs by location of firms
        psi_mat_II = m_psi_II * ones(n_grid_phi ^ 2, n_grid_del ^ 2); % average relationship among firms inland
        psi_mat_CC = psi_mat_II; % average relationship cost within the same region is assumed to be identical
        psi_mat_IC = m_psi_IC * ones(n_grid_phi ^ 2, n_grid_del ^ 2);
        psi_mat_CI = psi_mat_IC;
        psi_mat_F = m_psi_F * ones(n_grid_phi ^ 2, n_grid_del ^ 2);

        % Compute matrix of matching probabilities
        xi_max_S_inland_inland = max(pi_S_inland_inland ./ psi_mat_II, 0);
        xi_max_S_inland_coastal = max(pi_S_inland_coastal ./ psi_mat_IC, 0);
        xi_max_S_inland_Finland = max(pi_S_inland_Finland ./ psi_mat_F, 0);
        xi_max_S_inland_Fcoastal = max(pi_S_inland_Fcoastal ./ psi_mat_F, 0);
        xi_max_S_coastal_inland = max(pi_S_coastal_inland ./ psi_mat_CI, 0);
        xi_max_S_coastal_coastal = max(pi_S_coastal_coastal ./ psi_mat_CC, 0);
        xi_max_S_coastal_Finland = max(pi_S_coastal_Finland ./ psi_mat_F, 0);
        xi_max_S_coastal_Fcoastal = max(pi_S_coastal_Fcoastal ./ psi_mat_F, 0);
        xi_max_C_inland_inland = max(pi_C_inland_inland ./ psi_mat_II, 0);
        xi_max_C_inland_coastal = max(pi_C_inland_coastal ./ psi_mat_IC, 0);
        xi_max_C_inland_Finland = max(pi_C_inland_Finland ./ psi_mat_F, 0);
        xi_max_C_inland_Fcoastal = max(pi_C_inland_Fcoastal ./ psi_mat_F, 0);
        xi_max_C_coastal_inland = max(pi_C_coastal_inland ./ psi_mat_CI, 0);
        xi_max_C_coastal_coastal = max(pi_C_coastal_coastal ./ psi_mat_CC, 0);
        xi_max_C_coastal_Finland = max(pi_C_coastal_Finland ./ psi_mat_F, 0);
        xi_max_C_coastal_Fcoastal = max(pi_C_coastal_Fcoastal ./ psi_mat_F, 0);

        % [28]  MATCHING PROBS
        if matching_flag == "endogenous"
            m_S_inland_inland = G_xi(xi_max_S_inland_inland);
            m_S_inland_coastal = G_xi(xi_max_S_inland_coastal);
            m_S_inland_Finland = G_xi(xi_max_S_inland_Finland);
            m_S_inland_Fcoastal = G_xi(xi_max_S_inland_Fcoastal);
            m_S_coastal_inland = G_xi(xi_max_S_coastal_inland);
            m_S_coastal_coastal = G_xi(xi_max_S_coastal_coastal);
            m_S_coastal_Finland = G_xi(xi_max_S_coastal_Finland);
            m_S_coastal_Fcoastal = G_xi(xi_max_S_coastal_Fcoastal);

            m_C_inland_inland = G_xi(xi_max_C_inland_inland);
            m_C_inland_coastal = G_xi(xi_max_C_inland_coastal);
            m_C_inland_Finland = G_xi(xi_max_C_inland_Finland);
            m_C_inland_Fcoastal = G_xi(xi_max_C_inland_Fcoastal);
            m_C_coastal_inland = G_xi(xi_max_C_coastal_inland);
            m_C_coastal_coastal = G_xi(xi_max_C_coastal_coastal);
            m_C_coastal_Finland = G_xi(xi_max_C_coastal_Finland);
            m_C_coastal_Fcoastal = G_xi(xi_max_C_coastal_Fcoastal);
        else
            m_S_inland_inland = io_param_D; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_S_inland_inland(xi_max_S_inland_inland == 0) = 0;
            m_S_inland_coastal = io_param_D; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_S_inland_coastal(xi_max_S_inland_coastal == 0) = 0;
            m_S_inland_Finland = io_param_F; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_S_inland_Finland(xi_max_S_inland_Finland == 0) = 0;
            m_S_inland_Fcoastal = io_param_F; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_S_inland_Fcoastal(xi_max_S_inland_Fcoastal == 0) = 0;
            m_S_coastal_inland = io_param_D; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_S_coastal_inland(xi_max_S_coastal_inland == 0) = 0;
            m_S_coastal_coastal = io_param_D; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_S_coastal_coastal(xi_max_S_coastal_coastal == 0) = 0;
            m_S_coastal_Finland = io_param_F; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_S_coastal_Finland(xi_max_S_coastal_Finland == 0) = 0;
            m_S_coastal_Fcoastal = io_param_F; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_S_coastal_Fcoastal(xi_max_S_coastal_Fcoastal == 0) = 0;

            m_C_inland_inland = io_param_D; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_C_inland_inland(xi_max_C_inland_inland == 0) = 0;
            m_C_inland_coastal = io_param_D; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_C_inland_coastal(xi_max_C_inland_coastal == 0) = 0;
            m_C_inland_Finland = io_param_F; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_C_inland_Finland(xi_max_C_inland_Finland == 0) = 0;
            m_C_inland_Fcoastal = io_param_F; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_C_inland_Fcoastal(xi_max_C_inland_Fcoastal == 0) = 0;
            m_C_coastal_inland = io_param_D; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_C_coastal_inland(xi_max_C_coastal_inland == 0) = 0;
            m_C_coastal_coastal = io_param_D; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_C_coastal_coastal(xi_max_C_coastal_coastal == 0) = 0;
            m_C_coastal_Finland = io_param_F; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_C_coastal_Finland(xi_max_C_coastal_Finland == 0) = 0;
            m_C_coastal_Fcoastal = io_param_F; % .* ones(n_grid_phi^2, n_grid_del^2);
            m_C_coastal_Fcoastal(xi_max_C_coastal_Fcoastal == 0) = 0;
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%
        % Get the updated Phi and Del
        PHI_INLAND = kron(ones(n_grid_phi, n_grid_del), Phi_inland);
        PHI_COASTAL = kron(ones(n_grid_phi, n_grid_del), Phi_coastal);

        DEL_INLAND = kron(ones(n_grid_phi, n_grid_del), Del_inland);
        DEL_COASTAL = kron(ones(n_grid_phi, n_grid_del), Del_coastal);

        Phi_new_inland = (kron(phi, ones(1, n_grid_phi))) .^ (sig - 1) .* Kappa .^ (1 - sig) .* inland_status_iter .* entry_status_iter ... % network efficiency with empty network
            + (net_open) .* (alp / mu) ^ (sig - 1) * (newmatmean(PHI_INLAND .* m_S_inland_inland, g_chi, 0) ... % inland firms matched with all domestic inland supply firms
            + (tau_D) ^ (1 - sig) * newmatmean(PHI_COASTAL .* m_S_inland_coastal, g_chi, 0) ... % inland firms matched with all domestic coastal supply firms
            + (tau_B * tau_D) ^ (1 - sig) * newmatmean(PHI_COASTAL .* m_S_inland_Fcoastal, g_chi, 0) ... % inland firms matched with all foreign coastal supply firms
            + (tau_B * tau_D * tau_DF) ^ (1 - sig) * newmatmean(PHI_INLAND .* m_S_inland_Finland, g_chi, 0)); % inland firms matched with all foreign inland supply firms

        Del_new_inland = mu ^ (-sig) * (kron(ones(n_grid_del, 1), del')) .^ (sig - 1) .* inland_status_iter .* entry_status_iter ... % network demand with empty network
            + (net_open) .* mu ^ (-sig) * alp ^ (sig - 1) * (newmatmean(DEL_INLAND .* m_C_inland_inland, g_chi, 0) ... % inland firms matched with all domestic inland customer firms
            + (tau_D) ^ (1 - sig) * newmatmean(DEL_COASTAL .* m_C_inland_coastal, g_chi, 0) ... % inland firms matched with all domestic coastal customer firms
            + (tau_B * tau_D) ^ (1 - sig) * newmatmean(DEL_COASTAL .* m_C_inland_Fcoastal, g_chi, 0) ... % inland firms matched with all foreign coastal customer firms
            + (tau_B * tau_D * tau_DF) ^ (1 - sig) * newmatmean(DEL_INLAND .* m_C_inland_Finland, g_chi, 0)); % inland firms matched with all foreign inland customer firms

        Phi_new_coastal = (kron(phi, ones(1, n_grid_phi))) .^ (sig - 1) .* Kappa .^ (1 - sig) .* coastal_status_iter .* entry_status_iter ... % network efficiency with empty network
            + (net_open) .* (alp / mu) ^ (sig - 1) * (newmatmean(PHI_COASTAL .* m_S_coastal_coastal, g_chi, 0) ... % coastal firms matched with all domestic coastal supply firms
            + (tau_D) ^ (1 - sig) * newmatmean(PHI_INLAND .* m_S_coastal_inland, g_chi, 0) ... % coastal firms matched with all domestic inland supply firms
            + (tau_B * tau_DF) ^ (1 - sig) * newmatmean(PHI_INLAND .* m_S_coastal_Finland, g_chi, 0) ... % coastal firms matched with all foreign inland supply firms
            + (tau_B) ^ (1 - sig) * newmatmean(PHI_COASTAL .* m_S_coastal_Fcoastal, g_chi, 0)); % coastal firms matched with all foreign coastal supply firms

        Del_new_coastal = mu ^ (-sig) * (kron(ones(n_grid_del, 1), del')) .^ (sig - 1) .* coastal_status_iter .* entry_status_iter ... % network demand with empty network
            + (net_open) .* mu ^ (-sig) * alp ^ (sig - 1) * (newmatmean(DEL_COASTAL .* m_C_coastal_coastal, g_chi, 0) ... % coastal firms matched with all domestic coastal customer firms
            + (tau_D) ^ (1 - sig) * newmatmean(DEL_INLAND .* m_C_coastal_inland, g_chi, 0) ... % coastal firms matched with all domestic inland customer firms
            + (tau_B * tau_DF) ^ (1 - sig) * newmatmean(DEL_INLAND .* m_C_coastal_Finland, g_chi, 0) ... % coastal firms matched with all foreign inland customer firms
            + (tau_B) ^ (1 - sig) * newmatmean(DEL_COASTAL .* m_C_coastal_Fcoastal, g_chi, 0)); % coastal firms matched with all foreign coastal customer firms

        xi_bar_C_inland_inland = E_xi(xi_max_C_inland_inland);
        xi_bar_C_inland_coastal = E_xi(xi_max_C_inland_coastal);
        xi_bar_C_inland_Finland = E_xi(xi_max_C_inland_Finland);
        xi_bar_C_inland_Fcoastal = E_xi(xi_max_C_inland_Fcoastal);

        xi_bar_C_coastal_inland = E_xi(xi_max_C_coastal_inland);
        xi_bar_C_coastal_coastal = E_xi(xi_max_C_coastal_coastal);
        xi_bar_C_coastal_Finland = E_xi(xi_max_C_coastal_Finland);
        xi_bar_C_coastal_Fcoastal = E_xi(xi_max_C_coastal_Fcoastal);

        RC_inland = newmatmean(xi_bar_C_inland_inland .* psi_mat_II, g_chi, 0) + ...
            newmatmean(xi_bar_C_inland_coastal .* psi_mat_IC, g_chi, 0) + ...
            newmatmean(xi_bar_C_inland_Finland .* psi_mat_F, g_chi, 0) + ...
            newmatmean(xi_bar_C_inland_Fcoastal .* psi_mat_F, g_chi, 0);

        RC_coastal = newmatmean(xi_bar_C_coastal_coastal .* psi_mat_CC, g_chi, 0) + ...
            newmatmean(xi_bar_C_coastal_inland .* psi_mat_CI, g_chi, 0) + ...
            newmatmean(xi_bar_C_coastal_Finland .* psi_mat_F, g_chi, 0) + ...
            newmatmean(xi_bar_C_coastal_Fcoastal .* psi_mat_F, g_chi, 0);
        RC_total = RC_inland + RC_coastal;

        % Store outputs
        out.Phi_inland = Phi_inland;
        out.Del_inland = Del_inland;
        out.Phi_coastal = Phi_coastal;
        out.Del_coastal = Del_coastal;

        out.Phi_new_inland = Phi_new_inland;
        out.Del_new_inland = Del_new_inland;
        out.Phi_new_coastal = Phi_new_coastal;
        out.Del_new_coastal = Del_new_coastal;

        out.RC_inland = RC_inland;
        out.RC_coastal = RC_coastal;
        out.RC_total = RC_total;

        out.m_S_inland_inland = m_S_inland_inland;
        out.m_S_inland_coastal = m_S_inland_coastal;
        out.m_S_coastal_inland = m_S_coastal_inland;
        out.m_S_coastal_coastal = m_S_coastal_coastal;

        out.m_S_inland_Finland = m_S_inland_Finland;
        out.m_S_inland_Fcoastal = m_S_inland_Fcoastal;
        out.m_S_coastal_Finland = m_S_coastal_Finland;
        out.m_S_coastal_Fcoastal = m_S_coastal_Fcoastal;

        out.m_C_inland_inland = m_C_inland_inland;
        out.m_C_inland_coastal = m_C_inland_coastal;
        out.m_C_coastal_inland = m_C_coastal_inland;
        out.m_C_coastal_coastal = m_C_coastal_coastal;

        out.m_C_inland_Finland = m_C_inland_Finland;
        out.m_C_inland_Fcoastal = m_C_inland_Fcoastal;
        out.m_C_coastal_Finland = m_C_coastal_Finland;
        out.m_C_coastal_Fcoastal = m_C_coastal_Fcoastal;
    end % PhiDel function ends

    % [3] Definition of faIter function
    function out = faIter(fa_iter, Del_sol, Del_H_iter)
        fa_new = max((((t * eps * (sig - 1) * (kron(phi, ones(1, n_grid_del))) .^ (sig - 1) .* (fa_iter .^ (sig - 2)) * (mu - 1) * Del_H_iter .* Del_sol) .^ (1 / sig)) - t * eps) ./ w, 0.0001);
        out.fa_new = fa_new;
    end % End of fa function

    
    %[4] Definiation of matmean function
    function Dbar = matmean(D, g, renorm)
        % set default value of renorm
        if nargin == 2
            renorm = 0;
        end

        % renormalize if needed
        if renorm == 1
            g = g / sum(g(~isnan(D)));
        end

        % compute product of data and density
        Dg = D .* g;

        % compute mean
        Dbar = sum(Dg(~isnan(Dg)));

    end

    %[5] Calculate the weighted average of each block from kronecker products
    %(matrix of matmeans)
    function output = newmatmean(x, g, renorm)
        % set default value of renorm
        if nargin == 2
            renorm = 0;
        end

        % renormalize if needed
        if renorm == 1
            g = g / sum(g(~isnan(D)));
        end

        G = kron(ones(n_grid_phi, n_grid_del), g);
        G = x .* G;
        output = squeeze(sum(sum(reshape(G, n_grid_phi, n_grid_del, n_grid_phi, n_grid_del), 1), 3));

    end % End of newmatmean function

end % End of the entire function
